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We consider the classical response of a strongly chaotic Hamiltonian system. The spectrum of such 
a system consists of discrete complex Ruelle-Pollicott (RP) resonances which manifest themselves 
in the behavior of the correlation and response functions. We interpret the RP resonances as the 
cigenstates and eigenvalues of the Fokker-Planck operator obtained by adding an infinitesimal noise 
term to the first-order Liouville operator. We demonstrate how the deterministic expression for 
the linear response is reproduced in the limit of vanishing noise. For the second-order response we 
establish an equivalence of the spectral decomposition with infinitesimal noise and the long-time 
asymptotic expansion for the deterministic case. 
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I. INTRODUCTION 

Vibrational molecular motion is an interesting and im- 
portant example of complex dynamics. Due to nonlincar- 
ity, nonadiabatic effects and interactions with the envi- 
ronment, the behavior of such systems is remarkably rich, 
and their theoretical description away from the equilib- 
rium becomes a challenge. In many cases, e.g. when 
the frequencies of the vibrational modes do not exceed 
the temperature, a completely classical description is ad- 
equate. The dynamics far from equilibrium is typically 
irregular and can exhibit chaotic features. In the ex- 
treme case of strong chaos all trajectories in phase space 
are unstable, and the phase-space dynamics shows mix- 
ing features. 

Spectroscopic methods represent powerful tools for ob- 
taining detailed information on vibrational dynamics. 
The outcome of spectroscopic experiments can be nat- 
urally interpreted in terms of optical response functions, 
either in the time or in frequency domain. The multi-time 
(nonlinear) response functions that constitute the main 
objects in multi-dimensional time-domain spectroscopy 
are usually analyzed using their spectral decompositions, 
i.e. the Fourier transforms with respect to the time in- 
tervals. These time intervals can be viewed as the time 
delays between the exciting pulses, provided the latter 
are short compared to the typical time scales of the 
system dynamics. Periodic features of response func- 
tions are routinely interpreted as signatures of periodic 
motions. The central frequencies of diagonal and off- 
diagonal (cross) peaks are associated with some coherent 
vibrations, whereas the peak shapes contain information 
on the system-bath coupling. Smooth behavior of re- 
sponse functions in the case of almost harmonic vibration 
of the primary system is achieved due to its coupling to a 
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macroscopically large number of the bath motions. The 
bath is often considered in the harmonic approximation: 
such model is known in the spectroscopic literature as a 
multi-mode Brownian oscillator. 

The situation is totally different in the case of chaotic 
motion of the primary system. The spectrum of a 
strongly chaotic system is known to consist of so- 
called Ruelle-Pollicott (RP) resonances that represent 
the eigenmodes of the Perron- Frobenius operator [l|, 0] • 
The imaginary part of any particular resonance describes 
an oscillating feature in the system correlations, whereas 
its real part is responsible for the correlations decay. An 
individual RP resonance is not directly related to any 
particular periodic motion, although their positions can 
be expressed in terms of all periodic orbits in a very 
collective way via the dynamical ^-function. These col- 
lective resonances in chaotic systems should be distin- 
guished from the signatures of stable periodic motions. 

Nonlinear response functions of quantum systems can 
be conveniently represented in terms of spectral decom- 
positions via the system's stationary states. These spec- 
tral decompositions that are often referred to as Bloem- 
bergen's expressions allow to relate resonances in the 
spectroscopic measurement data to the transitions be- 
tween the stationary states. A natural question that 
arises in the context of approaching the classical limit 
is what would be the classical counterpart of the quan- 
tum Liouville space spectral decomposition for the linear 
and nonlinear response functions? In quantum mechan- 
ics the system state can be described by a density matrix, 
whereas the evolution is governed by the quantum Liou- 
ville operator. One can try to reach the classical limit 
by replacing the quantum density matrix by the classical 
phase-space distribution using the Wigner transforma- 
tion. The quantum Liouville operator should be natu- 
rally replaced by its classical counterpart. 

This conceptually straightforward approach, however, 
faces certain major difficulties. Quantum Liouville opera- 
tors are second-order elliptic operators, whose eigenstates 
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belong to certain natural Hilbert spaces. In particular in 
the case of compact (restricted) coordinate spaces the 
corresponding spectra are discrete, and finding spectral 
decompositions does not face conceptual difficulties, at 
least on the physical level of rigor. The situation with 
the classical limit is much more involved since the classi- 
cal Liouville operator is a first-order non-elliptic operator 
that describes the propagation of phase-space distribu- 
tions along classical trajectories. In particular any non- 
periodic classical trajectory generates a family of eigen- 
states of the classical Liouville operator, concentrated on 
the trajectory with all possible eigenvalues. Such a spec- 
trum contains no useful information on the system relax- 
ation. 

A meaningful spectrum that provides detailed infor- 
mation on the relaxation of a strongly chaotic system is 
known to be represented by the RP resonances. One of 
the ways to reproduce the RP resonances as eigenvalues 
of the classical Liouville operator is based on the ap- 
propriate simultaneous definition of the functional space 
where the operator acts. These functional Hilbert spaces, 
often referred to as rigged spaces, are very different from 
"standard" Hilbert spaces involved in spectral decompo- 
sitions of quantum operators. They are usually chosen 
on the case-to-case basis. 

To avoid these difficulties we follow a more physical ap- 
proach introduced in Ref. H. The approach is based on 
introducing weak Langevin noise, followed by considering 
the limit of vanishing noise. It has several major advan- 
tages. First of all, this describes a realistic situation since 
any system is at least weakly coupled to some environ- 
ment, and in many cases the system-bath coupling can 
be described on the level of Langevin noise. Second, it 
allows to avoid dealing with a problem of choice of the 
appropriate Hilbert space. The stochastic Langevin pro- 
cesses can be described by adding a diffusion operator 
to the Liouville operator L, which results in the Fokker- 
Planck operator C = -(k/2)V 2 + L. The Fokker-Planck 
operator £ is a second-order elliptic operator, and its 
spectral decomposition is free of the aforementioned dif- 
ficulties that arise in the case of classical Liouville opera- 
tor. The RP resonances are obtained from the spectrum 
of the Fokker-Planck operator by taking the limit k — > 0. 
This allows to interpret introducing infinitesimal noise as 
a regularization procedure similar to coarse graining. 

The dependence of the eigenfunctions of the Fokker- 
Planck operator £ on the noise strength n is nonanalyt- 
ical. In the limit k — » of the vanishing noise strength 
the smooth eigenfunctions turn into generalized functions 
(distributions). The rigged Hilbert spaces are reproduced 
from the "standard" Hilbert spaces in the limit k — > 0. 
They are spanned on the generalized eigenfunction ob- 
tained from the eigenfunctions of £(k) by applying the 
noiseless limit. 

The calculation of the linear response function in terms 
of the generalized functions does not pose a problem, 
since the expression involves an expansion of the smooth 
function over generalized functions. The calculation of 



nonlinear response functions is essentially more compli- 
cated because one needs to define expansions of the gen- 
eralized functions over similar generalized functions. We 
avoid dealing with the latter problem by performing all 
calculations for weak, yet finite noise, followed by apply- 
ing the limit k — > to the final expressions. 

In our previous work [J] we demonstrated the conver- 
gence of the second-order response function for strongly 
chaotic systems. As an example of the strong chaos, we 
considered free motion on a compact surface of constant 
negative curvature. The model has been studied since 
more than a century and served as a prototype for quan- 
tum chaos [E B ■ The motion on a surface of genus 
2 is in particular known as a Hadamard billiard. Using 
the dynamical symmetry (DS) of the system we found 
analytical expressions for the response functions. The 
long-time asymptotic expansion of the second-order re- 
sponse function turned out to have the form of the double 
decomposition in the resonances that appear in the linear 
response. 

In the present work we find the decomposition of the 
linear and nonlinear response in Rucllc-Pollicott reso- 
nances for this classical strongly chaotic model by intro- 
ducing weak Langevin noise. The spectral decomposition 
in this case is conceptually straightforward. We further 
demonstrate that the decomposition coefficients converge 
in the noiseless limit K — ► 0, and the resulting spectral 
decompositions reproduce the asymptotic expansions for 
the purely classical response functions computed in our 
earlier work Q. To the best of our knowledge, this is the 
first calculation of the classical nonlinear response for a 
chaotic flow that uses spectral decomposition based on 
the noise regularization. 

Our paper is organized as follows. In Sec. [Til we sum- 
marize the geometry and dynamical symmetry of free 
motion on a compact surface of constant negative curva- 
ture. Statistical description of the response is presented 
in Section IIIIl We regularize the dynamics by adding 
noise to the Liouville operator and derive general forms 
of the spectral decompositions in Section [TV] The eigen- 
states of the resulting Fokker-Planck operator are found 
in Section [Vj Explicit forms of spectral decompositions 
are obtained in Section I VII for the linear response and in 
Section IVlII for the second-order response. 



II. CHAOTIC MODEL SYSTEM: 
LIOUVILLE-SPACE PICTURE AND 
DYNAMICAL SYMMETRY 

Following 0| we consider free motion on a 2D compact 
surface M 2 of constant negative curvature (Gaussian cur- 
vature) . This strongly chaotic system is described by the 
classical free-particle Hamiltonian 



2m 



9 ik PzPk 



2m ' 



(1) 



3 



that depends on the absolute value £ of the momentum 
p only. The curvature of the 2D configuration space is 
expressed in terms of the metric tensor g lk . The Hamil- 
tonian classical dynamics preserves the smooth compact 
3D manifold M 3 that corresponds to a fixed value of en- 
ergy. Points x € M 3 of the reduced phase space are de- 
scribed by two coordinates r € M 2 and the momentum 
direction angle 6. Hereafter, we will use dimensionless 
units so that the mass m = 1 and the curvature K = — 1. 
Despite the complexity of the flow due to its chaotic na- 
ture, its strong dynamical symmetry (DS) enables an an- 
alytical treatment [!]• 

The geodesic flow f] = Lr), where 77 = (r,p) = (x,() 
denotes a point in phase space, is generated by the Liou- 
ville operator 



dH ck dH ck 
dp dr dr dp 
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Hereafter, we identify any vector field with the corre- 
sponding first-order operator of the derivative in the vec- 
tor field direction. We denote by a\ the vector field that 
determines the phase space velocity: L = £a\. We fur- 
ther introduce the second vector field a z = dj 88 in the 
tangent space, and finally set 02 = [cri, cr z ] (see e.g. Refs. 
13, @ for the details). A simple local calculation yields 
[<72) fz] = — <J\ arL d [<Ji,<72] = —Ka z , which implies that 
in the constant negative curvature case the vector fields 
1 , (T2 and a z form the Lie algebra so(2, 1). The group 
50(2, 1) action in the reduced phase space M 3 is ob- 
tained by integrating the so(2, 1) algebra action. DS 
with respect to the action of the group G = 50(2, 1) 
does not mean symmetry in a usual sense, i.e. that the 
system dynamics commutes with the group action, but 
rather reflects the fact that the vector field L that de- 
termines the classical dynamics is represented by an ele- 
ment of the corresponding Lie algebra so(2, 1), whereas 
the stable and unstable directions of our hyperbolic flow 
are determined by a z =p a%. The algebra generators 07, 
l,2,z are anti-Hermitian, i.e. satisfy the equalities 
= —01 with respect to the natural scalar product 
(ip,ip) = J dxip* (x; s)ip(x; s). In addition, the Poisson 
bracket of two functions f(x,() and g(x,() reads 
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The smooth action of G = 50(2,1) in M 3 can be 
interpreted as that the space H. of smooth functions in 
M 3 constitutes a representation of G, which turns out to 
be a unitary representation (see Refs. H, [^, [Toh . and 
therefore can be decomposed into a direct sum of ir- 
reducible representations of G. Stated differently, any 
distribution in the reduced phase space M 3 can be de- 
composed in irreducible representations. DS implies that 
the distributions in different representations evolve inde- 
pendently. We focus on the principal series represen- 
tations of 50(2, 1) since only these provide experimen- 



tally interesting contributions to the linear and second- 
order response [J]. A principal series representation H. s 
is labeled by an imaginary number s, Ims > 0. Eigen- 
functions ipk(x;s) of the momentum rotation operator 
<7 Z , hereafter also referred to as angular harmonics, form 
a convenient basis set in the irreducible representation 
TL S . The functions V'o( a; ;s) do not depend on the mo- 
mentum direction and can be viewed as eigenfunctions 
of the Laplacian operator on our compact Riemann sur- 
face M 2 [J, [l(| [H| • The Laplacian eigenvalues provide a 
set of numbers s € Spec(Af 2 ) according to the equation 
S7 2 ipo(x; s) — (s 2 — l/4)?/'o( a; ; s). The angular harmonics 
have the following properties: 



a z ip k {x;s) = ikt/j k (x;s) 



1 



(4) 



a±ij)k{x; s) = ( ±k + - - s ) ip k ±i(x; s) , 



where we introduced the raising and lowering operators 
a± = <7i ± ioi- The operators a± are anti-Hermitian 
conjugated: <j\_ = — cr_. 

The description of unitary representations of G @, Q 
allows us to identify the functions ipk(x;s) in reduced 
phase space with the corresponding functions on the cir- 
cle ^>k{u) for any s £ Spec(Af 2 ). In the case of the 
principal series Res = 0, Ims > 0, and ^(it) form an 
orthogonal normalized basis set with the natural scalar 
product. Therefore, the normalized functions ipk(x;s) 
are naturally associated with ^(u) = exp(ifcu). On the 
circle the algebra generators are represented by (see Q) 



a z = — , cri = smii- 
au au 

d 1 - 2s 



<7 2 



COS It 



du 



1 - 2s 
+ — 

sin it . 

1 - 2s 



■ cos u , 



(5) 
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a± = exp(±iw) ( ^i — 
du 



Using the dynamical symmetry (DS) of the problem, 
one can calculate the response functions. Two-point cor- 
relations that are related to the linear response via the 
fluctuation-dissipation theorem have been also consid- 
ered in Ref. [ll|. Detailed calculations of the response 
functions are presented in our previous paper Q. Since 
the zero harmonics ipo(x; s) are momentum-independent, 
the expansion of the dipole f(r) in irreducible represen- 
tations has the form 



/= B s ip (x;s) 

seSpcc (M 2 ) 



(6) 



where only the principle series contributions are retained. 
Typically the dipole is represented by a smooth distribu- 
tion, and, therefore, only a small number Nf of repre- 
sentations are relevant in the expansion. For the sake 
of the presentation clarity we focus on the case Nf = 1. 
Generalizations to Nf > 1 are straightforward and can 
be easily performed. 
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III. LANGEVIN STOCHASTIC DYNAMICS: 
FROM LIOUVILLE TO FOKKER-PLANCK 
EQUATION 

A. Liouville equation and response functions 

Strongly chaotic systems are characterized by compli- 
cated irregular dynamics [l2| . Exponential divergence of 
initially close trajectories and mixing make the descrip- 
tion in terms of individual trajectories complicated and 
inadequate. An alternative picture, based on Liouville 
representation of classical mechanics, involves propagat- 
ing phase space densities. Suppose that the system ini- 
tially occupies a certain region in the phase space which 
correspond to a range of initial phase space variables. 
In the course of the evolution of a chaotic system, the 
shape of the region becomes stretched along the unstable 
directions and contracted along the stable ones. 

The response functions characterize the system re- 
sponse to a time-dependent external driving field £(t) 
and constitute the basic outcome of spectroscopic mea- 
surements. The external field is coupled to the system via 
the polarization /(r), also often referred to as the dipole. 
The function f(r) of the system coordinates is the clas- 
sical counterpart of the polarization operator. The total 
classical Hamiltonian of the system coupled to the driv- 
ing field has the form 

H T = H- f£(t) . (7) 

The initial equilibrium phase space distribution po 
starts to evolve once the external field is turned on. The 
phase space density p(t], t) at time t depends on the val- 
ues of the external field at preceding times. In most cases 
the observed signal at time t is directly related to the 
same polarization f(r), which describes the system-field 
coupling. Within the classical mechanics formalism the 
observed polarization is given by J drj f (r) p(rj , t) . The 
evolution of the phase space distribution p is governed 
by the classical Liouville equation 

(d t + L)p = £(t){f,p}, (8) 

where the Liouville operator L acts as the Poisson bracket 
with the Hamiltonian according to Eq. ([2|) . 

The sign of the Poisson bracket is defined by a conven- 
tion that {pi,rj} — Sij for canonically conjugate coordi- 
nates r and momenta p. Note that the Liouville operator 
is a first-order linear differential operator. As stated ear- 
lier, it can be identified with the vector field L = in 
the phase space. The Liouville operator is anti-Hermitian 
with respect to the conventional scalar product, so that 
L* = —L. This makes the evolution operator e~ Lt uni- 
tary and leads to the conservation of phase space volume 
(Liouville theorem). 

One can solve Eq. ([8]) iteratively considering the right- 
hand side as a small perturbation and taking p for the 
zero approximation. Then the phase space density, and 
hence the response (observed signal) can be naturally 



represented in the form of the functional Taylor expan- 
sion over the external field. The n-th order response 
function S n is defined as the coefficient (more precisely, 
an integral kernel) in front of the n-th power of the driv- 
ing field. The response function S n depends on n + 1 
time moments, and if the unperturbed system is initially 
at equilibrium, it actually depends on n time intervals. 

For the linear and second-order response functions we 
have the following expressions: 

S«(fi) = J dr,f( V )e- Lt Hf(v),Po}, (9) 

S (2) {hM) = J dvf(v)e- Lt2 {f(v),e- Lt Hf(v),Po}}- 

(10) 

Here e _Lt is the the phase space density evolution op- 
erator, often referred to as the Perron-Frobenius opera- 
tor. We understand e~ Lt g{rf) as the solution p(r), t) of 
the equation (dt + L)p = with the initial condition 
p(x,0) = g(x). Since L is a first-order differential opera- 
tor, the solution can be readily found using the method of 
characteristics. In other words, the first-order operator 
L allows to write e~ Lt g{rf) = g{eT Lt r]). 

B. Langevin dynamics and Fokker-Plank equation 

Phase space trajectories found from the Hamilton 
equations of motion are invariant with respect to the time 
reversal. This property is easily observed in the behavior 
of integrable systems. Namely, any integrable dynam- 
ics can be represented by a set of quasiperiodic motions, 
which implies reversibility and recurrence in the values of 
physical observables. The behavior of chaotic systems is 
quite different. Although described by the same formal- 
ism based on Eqs. ([2|) and ((8]) that possess time-reversal 
symmetry, one observes obviously irreversible features 
such as relaxation phenomena. 

The apparent paradox and its solution are well known 
in statistical mechanics. In the course of evolution of a 
chaotic system, more and more fine features develop in 
the phase space distribution. The distribution width de- 
creases exponentially along the stable directions. At the 
same time, physical quantities are represented by smooth 
functions in phase space. Therefore, fine features of the 
distribution are actually irrelevant for the smooth ob- 
servables. In particular, this results in the exponential 
damping of the linear and nonlinear response functions 
in a strongly chaotic system Q ■ 

Only coarse grained properties of the phase space den- 
sity remain relevant. For instance, in the long-time 
asymptotic, the strongly chaotic system may be found 
in the neighborhood of any point of the phase space with 
equal probability, which is reflected in the homogeneity 
of the stationary phase space density. A physically use- 
ful and meaningful definition of the evolution operator 
which relates initial and final distributions must rely on 
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some kind of regularization that eliminates unnecessary 
distracting details. Coarse graining can be introduced 
either directly in phase space or, alternatively, in the 
functional space of phase space distributions. Sometimes 
coarse graining is inevitable, e.g. in computer simula- 
tions, where numerical errors determine the precision. 

The finest scale of classical dynamics is limited by the 
onset of quantum effects. One of well-known examples 
is the quasiclassical calculation of the entropy. Quan- 
tum effects were also shown to remove unphysical power- 
law divergences in the nonlinear response functions of 
integrable systems [l3|, [lij. It has been demonstrated 
that the limit h — ► should be taken after calculating 
the long-time asymptotics. This property resembles the 
characteristic feature of our approach to the second-order 
response in a chaotic system, as presented below. 

Besides, real physical systems are never isolated. The 
influence of an environment shows up as noise or random 
forces at the level of equations of motion, and as un- 
observable degrees of freedom and diffusive behavior of 
relevant variables in the statistical description. Although 
in many cases interaction with the environment is suffi- 
ciently weak and can be neglected, it is often utilized as 
a convenient and physically meaningful way to perform 
calculations. 

The effect of the diffusion in the chaotic system is to 
regularize the long-time dynamics and to introduce irre- 
versible leveling of gradients in the stable direction. We 
regularize the operator generating the density evolution 
by adding a small diffusive term in the form of the second- 
order differential operator: 

t = o x - \o\ , (11) 

Since a z generates rotations of the momentum, k/2 > 
has the meaning of the angular diffusion coefficient. 
Factor 1 /2 is chosen for reasons of convenience. We will 
be interested in the case n <C 1. As we will see, the limit 
k — > does exist in the spectral decomposition, therefore 
in the limit the factor does not play any role. 

The resulting Fokker-Planck operator C defines the 
evolution of the phase space density. This statistical 
description for the density is equivalent to the descrip- 
tion in terms of the Langevin equations for the dynam- 
ical variables. Such equations contain a random force 
which tends to change the momentum direction keeping 
the energy constant. This can be formally described by 
a stochastic Liouville operator 

L st (t) = C(a 1 + 1 (t)a z ) (12) 

where j(t) is a random Gaussian Markovian process with 
the zero mean. The noise intensity is determined by n, 
and the two-point correlation function reads 

(7(*)7(0>=«*(*-0- (13) 
The addition of the diffusion term with the second 
derivative defines the spectrum of the resulting opera- 
tor in the space of smooth functions. The eigenfunc- 
tions become regular, differentiable functions. As we will 



see, they are still singular at k = 0, and in the limit 
k — ► they turn into generalized functions, which is nat- 
ural for the case of a small parameter in front of the 
highest derivative. 

The mixing property of chaos makes the type of noise 
irrelevant. Random force exerted on any of the phase 
space variables affects other variables if they are not fixed 
by conservation laws. Thus, in the reduced phase space 
represented by the shell of constant energy the mixing 
leads to the fast and irregular randomization of the po- 
sition variables r. 

As a result, the Liouville equation ||SJ) that describes 
the phase space density evolution in a system perturbed 
by an external field £, is replaced by the Fokker-Planck 
equation: 

(d t + (£) P = £(t){f,p}, (14) 

The evolution operator in the unperturbed system can be 
written symbolically as e~^ ct , in a full analogy with the 
noiseless case. Similarly, the iterative solution of Eq. (fl"4"|) 
yields the phase space density in a form of an expansion 
in powers off. The response functions are obtained from 
Eqs. © by replacing L with (£: 

SM(t) = J d(afe- CCt f-Po), (15) 
S^(h,t 2 ) = J dCC(fe- CCt2 f-e-^f-p Q ) , (16) 

where angular brackets stand for the integral over the 
reduced phase space, and we defined the action of the 
operator /_ on a function g(rj) as the Poisson bracket of 
/ and g(r]) so that f_g = {f,g}. 

IV. SPECTRAL DECOMPOSITION OF 
RESPONSE FUNCTIONS: GENERAL 
FORMALISM 

In the well-known situation with quantum response, 
the response functions can be readily represented in the 
form of spectral decompositions, since the infinitesimal 
evolution is determined by a Hermitian operator. The 
conjugation property with respect to a simple scalar 
product facilitates expansions in the basis of the eigen- 
states. 

We have a natural scalar product for functions in M 3 
which allowed to implement the unitary presentations in 
3D space by functions on the circle. With respect to this 
scalar product L is an anti-Hermitian first-order differ- 
ential operator. However, the Fokker-Planck operator C 
is obtained by adding a second-order Hermitian contri- 
bution to L. This makes the resulting operator C neither 
Hermitian nor anti-Hermitian. Namely, its adjoint is 

£t = crj - \al . (17) 
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It seems that there is no natural choice of the scalar prod- 
uct that would make our Fokker-Planck operator simply 
related to its adjoint. 

The expressions for the response functions include in- 
tegrations over the reduced phase space. The convolu- 
tion (tp(x)ijj(x)) = J dx ip(x)ip(x) of two functions <p{x) 
and ip(x) can be recast in the form of the natural scalar 
product as ((pip) = ( < P*,i>)- 

As shown by the construction of eigenmodes in Sec- 
tion [Vl the operator C is diagonalizable in H s , and 
the Jordan-block structures that are possible in a gen- 
eral case do not appear in our spectral decompositions. 
Therefore, any function in Ti. s can be expanded in the 
eigenf unctions ip\(x;s) that obey the equation 



(18) 



The expansion coefficients are linear functionals of 
tp{x; s) and hence can be represented by scalar products 
of tp(x; s) with certain functions denoted by <p\(x; s): 

<pfa s ) = ^(fifa; s M x i s))<Pnfa s ) • ( 19 ) 

A" 

It follows from the expansion of <px(x;s) and C.(f\{x;s) 
that the functions <p\(x;s) are the eigenf unctions of the 
operator & which satisfy the following properties: 



tfip\{x]s) = \*ip\(x;s), 
(V>*a{x; s)<p\(x; s)) =S f _ l \. 



(20) 
(21) 



We focus our detailed treatment on the simplest case 
Nf = 1 of a single irreducible representation contribut- 
ing to the dipole moment f(r). This corresponds to a 
single contribution in the expansion ([6]), and below we 
imply / = ipa(X) s). A generalization to an arbitrary lin- 
ear combination of Nf > 1 such terms with different s 
is straightforward for both response functions under con- 
sideration [J]. 

We can successively apply the expansion procedure 
(|19p and obtain any response or correlation function as 
a spectral decomposition. Based on Eq. (fT"5f we get the 
following form for the linear response function: 



•J \ 



(22) 



For the calculation of the second-order response start- 
ing from Eq. (|16p . we in addition introduce a double ex- 
pansion in the basis of the angular harmonics that are 
defined by Eqs. (g]): 



sW(t u t 2 ) = I dccJ2J2 e ^ ct2x ( 23 ) 

(<pM «/_ j> m ) (r m e- xctl Vx) (^xf-Po) ■ 

The double expansion in the angular harmonics yields 
a geometric matrix element (ip^ipoip m } studied in our pre- 
vious work 0. Note that/_ in the middle angular brack- 
ets contains a derivative d/dC, acting on all functions of 



the momentum to the right of it. This complication oc- 
curs due to the absence of the integration over £ in the 
angular brackets which correspond to projections onto 
the basis set in H s . 



V. EIGENSTATES OF FOKKER-PLANK 
OPERATOR 

As stated earlier the Fokker-Planck operator (fTTj) is 
neither Hermitian nor anti-Hermitian with respect to any 
natural scalar product and, therefore, not necessarily di- 
agonalizable. In what follows, we show that it actually 
is: We apply the representation on the circle to find its 
eigenstates and demonstrate that they constitute a basis 
set. 

The eigenfunctions and eigenvalues of C can be found 
by solving the ID eigenvalue problem 



£$a = A$ A . 



(24) 



on a circle with the Fokker-Planck operator [see Eqs. ([5]) 
and HI])] 



2 du 2 du 



1 - 2s 



(25) 



The first derivative in the operator can be eliminated 
by redefining the functions 



$ a (m) = e c ™"£ A (w). 



(26) 



In terms of the functions £a(w), the eigenvalue problem 
(|24|) assumes the form of a stationary Schrodinger equa- 
tion 



with the effective Hamiltonian 
k d 2 



H 



sin 2 u 



2 du 2 



2 k 



— s cos u . 



(27) 



(28) 



The Hamiltonian describes a quantum particle in a 
complex-valued potential on the circle. The operator 
in not Hermitian because of the imaginary part of the 
potential. However, since the Hamiltonian does not con- 
tain first derivatives, we can define a symmetric (non- 
Hermitian) scalar product V(£, <fi) — J du £(u)(j)(u) so 
that the Hamiltonian is self-adjoint with respect to it: 
V(fi^,<j)) — V(^,Ti4>). Eigenfunctions of H correspond- 
ing to different eigenvalues are orthogonal and can be 
normalized: 



2tt 



(29) 



The Hamiltonian H also possesses certain symmetries 
that simplify the analysis. First, it is invariant with re- 
spect to the sign change of u, and therefore all its eigen- 
functions are either even or odd functions of u. As we 
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will see, only even eigenstates contribute to the spec- 
tral decompositions of response functions. Therefore, we 
present the detailed expressions only for even eigenstates 
on the half-circle < u < tt. 

Since the potential contains a nonzero imaginary part, 
its eigenvalues can be complex. The second symmetry 
involves the complex conjugation of the Hamiltonian: 
7i*(u) = TL{u + 7r). Consequently, the eigenfunctions 
corresponding to the complex conjugate eigenvalues are 
related by 



Mu) = (&(u + 7r))* 



(30) 



Detailed analysis of the Schrodinger equation in the 
weak noise case k <C 1 is presented in Appendix [XJ The 
small parameter k allows to solve the Schrodinger equa- 
tion using the WKB method. The imaginary part of 
the potential energy is small compared to its real part. 
This supports the use of such terms as the "minimum" of 
the potential and the "under the barrier" wave function. 
Since the potential minima become deeper for smaller k, 
the eigenfunctions are concentrated near u = or u = tt. 
These states do not mix together since s is imaginary. 

Due to the compact nature of the circle, the spectrum 
of the Fokker-Planck operator is discrete. Its real part is 
positive and unbounded from above. For such large en- 
ergies that Re A 3> the spectrum can be adequately 
approximated by A„ ~ nv 2 /2 (with the eigenfunctions 
£a„ ~ e lvu ) of a free particle on the circle. 

In the limit k — > even low-energy eigenstates of the 
first set are concentrated near u = and have energies 



A„ = v 



1 

S+ 2' 



(31) 



with nonnegative even numbers v {y = 0, 2,4, . . .). This 
discrete equidistant spectrum of the Fokker-Planck op- 
erator in the noiseless limit k — » is a quite fascinating 
property. Infinitesimal noise regularizes the Liouvillian 
dynamics to yield a physical spectrum in the space of 
smooth functions. 

The normalized eigenfunctions that correspond to the 
specified above eigenvalues are concentrated near u = 0: 



A , v = H-V^V!)- 1 



-1/2 



(32) 
(33) 



where H v are Hermite polynomials. The eigenfunctions 



a COS u/ K 



£a„ (w) are singular at u = in the limit 



k — ► 0, as they should form a basis set of distributions 
with large gradients in the stable direction. 

In the limit k — > the eigenvalues in Eq. (f3"Tj) are 
determined by the single potential minimum at u — 0. 
The form (|32[) of the eigenfunctions is appropriate in the 
region where the potential may be approximated by a 
harmonic well. Outside the region, at \u\ > 1, the values 
of the eigenfunction are exponentially small and do not 
influence its normalization. 



A more rigorous treatment in Appendix provides the 
WKB approximation for the even eigenfunctions (u) on 
the whole circle. Under the potential barrier, if u > k 1 / 2 
and (tt—u) > k 1 / 2 the approximate solution at the energy 
A„ = v — s + 1/2 is found to be 

£(u) = A 2 ^2 s -i (sin^^cos^) 6 " 5 , (34) 
A 2 ^ = A , u e-i2» +x K -%. (35) 
In the vicinity of u — tt, at n — u <C 1, we obtain: 



(36) 



where 



V\ = v — Is = A s , 

2 

A 1 = - J Bii2"" 1 7r" 1/2 r(-^i)cos 



TKV\ 



Bi = A 2 „e~*2 



(37) 
(38) 
(39) 



Even low-lying eigenstates of the second set concen- 
trated near u — n have eigenenergies A* = u+s+1/2, and 
their eigenfunctions can be obtained by using Eq. (|30p 
which follows from the symmetry of the Hamiltonian. 

Eqs. (|32"l) - (f3"7| provide a zero-order approximation for 
the eigenfunction of the Hamiltonian (|28p on the circle. 
In the limit k < 1 one can distinguish two types of cor- 
rections to the specified above eigenvalues and eigenfunc- 
tions of low-lying states. The first one is due to the omit- 
ted terms in the potential in the vicinity of its minima. 
Such corrections are important in the treatment of states 
with v > 0. Corrections of the second type originate from 
the tunneling through the potential barriers and there- 
fore are exponentially small ~ e -2 / K . 

Details of computing the eigenfunctions are presented 
in Appendix [5J It turns out that a straightforward cal- 
culation of the higher-mode contributions to the spectral 
decomposition requires knowledge of eigenfunctions with 
increasing accuracy in k. In the next two sections we 
present a more elegant approach for which the approxi- 
mation given by Eqs. (|3l?|) - (f3"T|) is sufficient. 



VI. 



SPECTRAL DECOMPOSITION FOR 
LINEAR RESPONSE 



Spectral decompositions of the response functions in 
the eigenmodes of the evolution operator are obtained in 
Eqs. (|2"!il |2"3")) in the general form. The expressions are 
not really useful for the calculation, since it would in- 
volve several three-dimensional integrations that involve 
the functions not explicitly known. Nevertheless, the cal- 
culation can be carried out by implementing the repre- 
sentation in terms of functions on the circle. 
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We start with the spectral decomposition of the linear 
response, and employ the correspondence between the 
functions (p in M 3 and <I>(it) on the circle, to obtain the 
following expansion: 



$(«) = ^{V\^x{u) 



(40) 



Using the representation <&\{u) — e~ cosu / K {; X (u) and the 
orthonormality of the functions t;\ (u) with respect to the 
symmetric scalar product (|29|) . we can specify a rule to 
find the projection of ip on ip x : 



($Z<p) = / du<S>{u)e^£ x (u). 



(41) 



This implies that (p x is implemented on the circle by the 
function 



(42) 



To proceed with the calculation of the response func- 
tions, we introduce a notation 



R x , n = {(p\^ n ) = / due im e^ix{u) (43) 



for the coefficient of the expansion of ipn(x; s) in <p\(x; s). 
The expansion of <p x (x; s) in angular harmonics ip n (x] s) 
contains coefficients 



. du 

(1> n <px) = I ^ 



(44) 



The symmetry property (|30[) allows to relate them to 

R\,n'- 



(-D n 



(45) 



Therefore, the convolution in the first angular brackets 
in Eq. ([22]) is equal to {R\*,a)* /2tt. The action of the /_ 
in the last angular brackets, represented by a Poisson 
bracket ([3]), creates the following convolution 



(<p* x (Tiipo) = 1 ^ I ^ cosue"°^ £ x (u) , (46) 



where we have used the representation of a\ on the circle 
given by Eq. ([5]). Bearing in mind that ^ x is the eigen- 
function of Hamiltonian (|28p and integrating by parts, 
we find the following relation: 



{(p* x aii> ) = XR Xfi . 



(47) 



We conclude that the linear response function can be 
recast in the form of the time derivative of the two-point 



correlation function, in agreement with the fluctuation- 
dissipation theorem (FDT): 



oo 

S (1) (t) = ^ [ d( Qx , e- 



CAt^O (48) 



where we introduced 



Q\,n 



(-1)' 
2tt 



R\fi{R\- ,n)* • 



(49) 



Following the approach developed in Ref. [4j for the 
purely deterministic situation, we introduce the matrix 
elements of the evolution operator between n-th and zero 
harmonics: 



An((t;s)= I dxiP* n (x;s)e- ct M x ;s) ■ (50) 



Spectral decomposition 
An((t;s) 



-CAt 



(51) 



appears in the response functions of the system with 
noise. 

As we show below, the noiseless limit k — * of the 
series A n (t; s) reproduces the matrix element A n (t; s) of 
the deterministic evolution operator, calculated in Ref. 
|4j. These series and their coefficients Qa,u play an impor- 
tant role in the calculation of the second-order response 
function. 

We are now in a position to proceed with an explicit 
calculation of the coefficients in the spectral decomposi- 
tion of the linear response function. We focus on the first 
set of modes with the energies A„ = v — s + 1/2, whose 
eigenf unctions are concentrated in the neighborhood of 
u = for k <C 1. The coefficients Q Xy o for the other set 
of modes, with the energies A* = v + s + 1/2, can be 
then easily found by employing the symmetry, described 
by Eq. A 



CX',0 



{Qx,oT 



which also reflects the fact that the response function is 
real. 

The first two terms in the spectral expansion originate 
from the lowest RP resonances with the energies Ao = 
1/2 — s and Aq = 1/2 + s. The eigenfunction with A = 
1/2 — s, concentrated in the vicinity of u = 0, is given 
by £a ( u ) = (7tk) -1 / 4 exp(— u 2 /2k). The first integral, 
R\ ,o> i s calculated without much effort. This is done 
by noticing that the main contribution to the integral 
comes from |tt| < y/H. Expanding cos it = 1 -u 2 /2 in the 
exponential results in a Gaussian integral: 



R\ ,o 



due™'"" £ Xo {u) = e»A 



0,0 i 



(52) 



9 



where the normalization factor A$ q is related to k via 
Eq. (33]). 

The situation, however, is more complicated for higher 
modes with v > 0. Due to the orthogonality of the eigen- 
functions £\ and £a„ with v > 0, the integral vanishes if 
we employ the lowest-order approximation for the wave- 
functions: 



due * (;\(u) k, A Q _ v e* j due K H u [ —j= ) = . 



and we arrive at 



(53) 



This actually means that more accurate expressions are 
needed to find the first nonvanishing term. Note that one 
should use better approximations for both the exponen- 
tial e costl / K and the wave function. The corrections to 
the latter are computed by applying the standard quan- 
tum mechanical perturbation theory to the Schrodinger 
equation. The procedure is feasible for few low modes 
only, since the calculation for the higher modes requires 
higher orders of the perturbation theory. 

However, it is possible to overcome this difficulty and 
find an explicit expression for R\ v .a for all v by using the 
following trick. We notice that the dominant contribu- 
tion to the integral 



du e^ (cos u - ir /2 6w («) = 2- 2 »{-iy /2 A^ u e 1/K . 

(54) 

does not involve higher-order perturbative calculations 
described above. The integral is dominated by the re- 
gion \u\ < \/k, where we can approximate (cosu— l)"/ 2 = 
(— l/2)"> 2 u v and expand it in the Hermite polynomials 
Hj(u/ 't/k) with j = 0, 2, . . . , v. All terms except for the 
last one vanish due to orthogonality of the Hermite poly- 
nomials. Further approximations lead to negligible cor- 
rections in the limit n — > 0. We proceed by introducing 



Pi,> 



du i 



: (cosw - l) j £\(u) 



(55) 



Our goal is to calculate the integral -Ra„,o = Po,\ v - This 
will be achieved by expressing it via the known inte- 
gral P u /2,\„ for nCl. Consider XPj.x- Making use of 
A£ A = TL^x, we integrate in Eq. (155)) by parts, and neglect 
higher-order terms in k. This results in the recurrence 
relation 



P 



J-S + 



A - 2 j 



(56) 



which being applied vjl times and followed by setting 
A = \ v yields: 



Po,\„ = 



W2)!r(^) 



l-2s\ ■ P "/2,A„ 



(57) 



(_-\ \v/2-r I v+1- 
p _ a-1 -2v \ L > 1 V 2 



W2)!r(^) 



(58) 



The second integral (ipo^x^) m the linear response is 
expressed via (R\*,o)* [see Eqs. (gU) and IjlSjl]. Making 
use of the complex conjugation symmetry specified in 
Eq. ([30]) this can be recast as 



7T 

(R K ,o)* = J due-^i 



(59) 



The obtained integral looks more complicated compared 
to R\„fl, even for the lowest resonance with v = 0. We 
first point out the problem and present its straightfor- 
ward solution for v — 0. Afterwards we apply a more el- 
egant method, similar to the one developed above, which 
provides the solution for all modes. 

Before starting we emphasize that, due to the eigen- 
function symmetry, the discussion can be limited to the 
interval < u < ir. Three different analytical approxi- 
mate forms for the eigenfunction near u — 0, u — it, and 
under the barrier are specified in Eqs. ([32]) - (|37|) . As 
opposed to the case of -Ra„,0j the integrand in (i?>,j,o)* 
does not experience exponential dependence on u, and 
the exponential part of its dependence on k in the con- 
tributions from all three regions amounts to a common 
factor exp(— 1/k). 

In the region yfk <C (tt — u) <C 1, where the approxi- 
mate forms, given by Eqs. (j34|) and (f36|) match, the so- 
lution behaves ~ (ir — u)~ u ~ 1+2s . This region, where the 
two solutions match, provides an essential contribution 
to the integral. We can break up the integral into two 
parts by introducing an intermediate point 7r — u\ that 
belongs to the matching region. Since for small n ^ 1 
one can always choose k 1 / 2 < iii < 1, the approximate 
forms provide good approximations to the left and to the 
right of the intermediate point, respectively. We intend 
to show that the sum of two integrals is independent of 

Ml. 

The contribution from the region nearby u = 7r has the 
form: 



due " £\(u) — 



(60) 



e«V« J dz(A x e 2 ' > H Ul (z) + B x H- Vl -i{iz)) = 
o 

R i r _ nvxl2 H- Vl (wi/Vg -H- Vl (-wi/Vk) 

— E>\e*-yKe : — r . 

Avi sin(7r^i/2) 

We chose u\ so that k 1 / 2 « hi « k 1 / 4 <C 1 and re- 
tained first two terms in the expansion of cos u. The 
integral is calculated by expressing e~ z H Vl (z) in terms 
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of H- Vl -\(iz) and iz), since all three func- 

tions are solutions of the same second-order differen- 
tial equation (see Ref. 111). We also employ the solu- 
tion symmetry and the recurrence relation -j^H- Ul [iz) — 
—2iuiH— Vl —i(iz) for Hermite functions. 

The integration with the under-the-barrier function 
can be safely extended to u = 0, since the deviation of 
the approximate integrand from the actual function in a 
short region of length ~ y/n is finite. Then the integra- 
tion is performed exactly in terms of the hypergeometric 
function: 



7T — Hi 

du 

o 



/ . u \ v ( u\ ! 
{ 8m 2) { C0S 2) 



i_ r(^i)r(^; 
r (s + *) 



(61) 



2 2 F 1 (l,s+j^, S in 2 ^) 



v+1 



We further expand Eq. ^ in <Jk/ui < 1 and Eq. |[6"T]) 
in ui -C 1. In the sum of two pieces of the integral ([59]) . 
the strongest dependence on cancels out. For the low- 
est resonance with v — there are no other contributions 
comparable with the constant in Eq. (|61|) . which is not 
the case if v > 0. 

As v increases the divergence of the under-barrier solu- 
tion becomes stronger, and one should take into account 
all Mi-dependent terms with non-positive real parts of the 
exponents that appear in Eqs. (|60|) and (j6"Tj) . For large 
v this requires higher order of the perturbative calcula- 
tions. One might have argued that, since the corrections 
to the wave functions have the same structure as the ini- 
tial approximations for them, all relevant i^-dependent 
terms should cancel out, and the corrections to the con- 
stant are small in k. 

However, we can avoid exact unpleasant calculations 
by applying the method developed above. This is 
achieved by relating (R\* j0 )* to the integral that can 
be calculated using the principal approximation for the 
eigenfunction specified in Eqs. (|32l - I57j) . According to 
the definition ([55]) of Pj,x and due to the wavefunction 
symmetry §5U§, we have (ii^j.o)* = (-Po.Aj)*- The recur- 
rence relation (p)6")) yields 



(Po,> 



^)r( s -|) 



(_ 2 )i+-/2 r ( s + i)r(i + s) 



(Pi 



l+;v/2,A; 



(62) 



Next, we consider the integral 

(JWaaj)* = ( 63 ) 

{-l) 1+v / 2 / due-^(l + cosu) 1+l//2 a>) 



and notice that for its calculation it is sufficient to take 
the under-barrier approximation (|34[) for £\ u on the 



whole half-circle. The resulting integral can be easily 
calculated: 



A 2 ^(_X)l+^/22(^+2s+3)/2 r (~27 r ( 1 + S ) 



(64) 



r((2s + ^ + 3)/2) 



The integral converges, and the differences between the 
approximate and exact integrands in the narrow regions 
of width ~ \/k near u — and u — ir are smaller than 
~ A2mK u / 2 and ~ Ai^k sJv1 I 2 , respectively. The calcu- 
lation of the second integral is completed by combining 
Eqs. (|6"2"j) and (f6"4")) followed by expressing A 2 ^ v via Ao, v , 
according to Eq. (f35|) : 



(65) 



r ( v+1 \ v ( 2s ~ v ^ 
(R x . y = ,Ve-^~" /2 2 1+ ^- 1 2 } { 2 ' 



r la- 



The method is naturally also applicable to the case v = 0, 
which has been treated earlier, using a more straightfor- 
ward approach, with much more effort involved. 

We complete the calculation by inserting Eqs. ([55]) and 
([65]) into Eq. (j49|) to find that Qa„,o and Qaj.o have finite 
limits at n — > 0, namely 



lim Q A o 
lim Q A *,o 



= tt(i//2)! r(i + 
(hmQA^o)*. 



COt(7Ts) 



(66) 



We conclude the section by demonstrating that the spec- 
tral decomposition of the linear response in the noiseless 
limit k — > reproduces the asymptotic expansion of its 
purely deterministic counterpart presented in Ref. 0. To 
that end we show that in the noiseless limit <2a„,o and 
Qa* ,o coincide with the coefficients in the expansion 



A n (t;s)= ( Q»,ne st + Qv.ne~ st ) e 

j/=0,2,... 



st) _-(i/+l/2)i 



(67) 



at n = 0, where A n (t;s) is defined as a matrix element 
of the Liouvillian deterministic evolution: 



A n {(t;s)= I dx%b* n (x;s)e Lt ^ Q {x; s) 



(68) 



The integral can be calculated using the representa- 
tion on the circle by solving the ID Liouville equation 
dte-^Wolu) = -a^-^oiu) [see Refs. HQ: 



, , \ 2(e- 2 *-l)"r(n + i- S ) , /2 
A n (t;s) = -i J_„ ; 2 V*/ 2 x 



Re> 



%^r(i - s ) 

r(s)e s * „ / . 1 1 



(69) 



— — t— - 2 F 1 [n+--s,n+-,l-s,e 24 
T(n+i + s) V 2 2 



The expression is substantially simplified in the relevant 
for the linear response case n = 0. Expanding the Gauss 
hypergeometric function 2 F\ in A (t; s) into the hyper- 
geometric series of e~ 2t to obtain Q v $ and Q^o, we see 
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directly that they are indeed reproduced by the noiseless 
limit of the coefficients Qa„.o an d Qa;,o given by Eq. (f6"6"j) . 
and in the limit k — > 0, i.e. 



Q 



and 



Qa*0 — Qufi ■ 



(70) 



VII. SPECTRAL DECOMPOSITION FOR THE 
SECOND-ORDER RESPONSE: NOISE 
REGULARIZATION 

Spectral decomposition of nonlinear response functions 
in the case of finite noise is conceptually straightforward. 
The general expression for S^{t\, £2), decomposed in 
the eigenmodes of the Fokker-Planck operator C(k) is 
given by Eq. ([2"3"]) . In this section we demonstrate that 
in the noiseless limit k — > the expansion coefficients 
converge to the coefficients of the long-time asymptotic 
expansion of the purely deterministic k = response that 
have been derived in our previous work For the sake 
of simplicity we focus on the Nf — 1 case when only 
one irreducible representation contributes to the dipole, 
i.e. / = ipo(x;s). The expression ([2"3"]) for the second- 
order response contains four matrix elements calculated 
in the previous section. The fifth matrix element in 
the middle angular brackets requires a careful treatment 
since /_ includes the operator a-^od/dC, that acts on all 
momentum-dependent functions to the right. We first 
perform integration over the reduced phase space, repre- 
sented by the middle angular brackets, which results in 
a C-dependent expression that also includes derivatives. 
Integrating over £ by parts, and employing the symme- 
tries, we obtain the second-order response function in the 
form of the following spectral decomposition: 



n=0 



(",- - a,, hi) ( ( n + s + - 



where Qa,m is given by Eq. 
a_„ coefficients 



(fit 2 + n 



(71) 



The symmetric a n = 



<i„~ I dx %p* n (x; s)ip (x; s)ip n (x; s) (72) 

' M 2 

are purely geometrical factors that do not depend on the 
dynamics and are Q ■ They are all related to eto with the 
help of the recurrence relations 



8n 2 + 1 - 4s 2 (2n - l) 2 - 4s 2 



(2n+l) 2 -4s 2 (2n+l) 2 -4s : 



O n -l (73) 



and, therefore, can be expressed in terms of the Laplacian 
operator eigenfunctions ipo(x; s). 



The spectral decomposition [7T] is an exact expression 
for any given k. The diffusion coefficient K enters the ex- 
pression via the eigenvalues as well as via <2a,ti, expressed 
in terms of the eigenfunctions as specified in Eq. ([ST]) . In 
what follows we determine the leading contributions in 
the noiseless limit k — + 0. 

Compared to the linear case, noise plays a more del- 
icate role for the spectral decomposition of the second- 
order response function. It provides convergence of the 
series over angular harmonics in the expression for the 
response function that does not appear in the linear case 
[compare Eqs. (gSJ) and ([71])] ■ 

To retrieve the asymptotic behavior of Qa,™ for k«1 
and n>lwe apply the method developed in Section [VTl 
to derive the following recurrence relation: 



A — ) />\.„ - 



(74) 



2n + 1 - 2s 



-R\ 



2n - 1 + 2s 



71+1 



R\,n-1 



Specifically we replace A£a with Tl^\ in the integral rep- 
resentation [Eq. ([4*3"]) ] for \R\, n followed by integrating 
by parts. Eq. (|74p results in another exact recurrent re- 
lation 



Qx,n 



+ 1 



2n - 1 - 2s „ 4A - 2^n 2 ^ . , 

;Qa,«-i~ n _ , ; , n Qx,n- (75) 



2n + 1 + 2s ' 



2n + 1 + 2s ' 



Since only symmetric eigenfunctions contribute to the 
spectral decomposition, we have Q\,- n — Q\, n , which 
allows to express all quantities Qa.k via <2a,o specified 
by Eqs. @gj>. 

We start with the limit k —> for fixed n. Although 
the dependence of R\^ n on k is singular, the coefficients 
Qx.n have well-defined limits at k — > that are equal 
to the corresponding coefficients in the expansion (|67p of 
the purely deterministic matrix elements A n (t; s). 

This can be established in the following way. Viewing 
A n (t; s) as a matrix element of the evolution operator 
between zero and n-th harmonic, we employ the identity 
(j 1 = (<r + +<7_) and the relations between the neighboring 
harmonics that follow from Eqs. ([5]) to express dtA n {t; s) 
in terms of A n -i(t;s) and A n+ i(t;s). Finally, the re- 
sult is expanded in e~ 2t according to Eq. ([67]) . For the 
components oscillating as oc e s */ 2 we arrive at a relation: 



_2n+l + 2s 2n-l-2s 



'u,n—l 1 

(76) 



which that can be viewed as the noiseless limit of the 
recurrence relation (|75p . Combined with the already es- 
tablished equivalence lim K _^o Qa„,o = Qv,o f° r the- zero 
term, this implies the equivalence lim K ^o Q\ v .n = Qv,n 
for any given n. 

So far the equivalence has been established for the RP 
resonances with \ w = v — s + 1/2. The coefficients for 
the other set of RP resonances can be easily found by 
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employing the symmetry of the recurrence relation (175 
combined with Qa*.o — (Qa o)*> which results in: 



r(n + i + «)r(i - s) 



(77) 



in the limit k — ► 0. This establishes the equivalence 
lim K ^ 2a*, n = <9i/,n f° r the other set of resonances. 

We are now in a position to demonstrate that the series 
in angular harmonics n for any coefficient in the spectral 
decomposition converges for small < n <C 1 and in the 
noiseless n — > limit reproduces the corresponding coef- 
ficient in the long-time asymptotic series for the purely 
deterministic response function S^(t 1 ,t 2 ). 

Explicit expressions for the coefficients Q v ,n and Qv,n 
that enter the expansion (|67[) for A n (t; s) become increas- 
ingly lengthy as v grows. The coefficients can be obtained 
by expanding the deterministic evolution matrix element 
in Eq. (|69[) in powers of e~ 2t . The simplest expressions 
can be found for the lowest modes with v = whose 
energies are A = ±s + 1/2. For example, for the mode 
with A = Aq=s + 1/2, the recurrence relation (|75|) im- 
plies Qo,n = (-l) n Qo,n, where Qo,n = Qq u and Q , n is 
specified in the r.h.s. of Eq. (JBSJ) taken at v = 0. 

To analyze the limit n — > oo for fixed finite n <^ I 
we view Rx,n as the Fourier coefficients of the smooth 
function $\(u). Consequently, Q\ :7l decay faster than 
any power of n for n — > oo, and indeed we find from the 
recurrence relations the intermediate asymptotic of Q\. n 
at 1< n < k^ 1 : 



Qx,n oc (-1) 1 



2 /4 



(78) 



For larger n the recurrence relation implies 
QA,n/QA,n-i = V( Kn )j which also leads to the de- 
cay faster than any power law as n increases. Finally we 
compare Eq. (f78|) with the asymptotic n> 1 form of 
the deterministic coefficients Q v n to derive 
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(79) 



which is valid in the range 1 <C n <C k^ 1 for small k. 
Eq. (j?9")) means that noise provides a homogeneous cut-off 
for all RP resonances and the series over angular harmon- 
ics for the spectral decomposition coefficients converges 
atn~ kT 1 ! 2 < xT 1 where Eq. (J75J) holds. 

It remains to be demonstrated that the resulting spec- 
tral decomposition in the n —> limit of vanishing noise 
coincides with the asymptotic long-time expansion of 
S( 2 '(tx,t2) in the deterministic case. The latter expan- 
sion has been derived in Ref. 0. 

First of all, we observe that the eigenvalues of the 
Perron-Frobenius operator converge in the k — > to the 
factors that appear in the expansion of A n (t; s) . The 
same factors that appear in the spectral decomposition 
of linear response have been interpreted in Ref. [ll[ as 
RP resonances. While the linear response function may 
be represented in the form of a converging series, the ex- 
pansion of appears to be more involved. If k = 0, 



for a pair of resonances, the series (fTTj) over n diverges, 
as clearly seen from the power-law growth in Eq. (|78D . 

In the deterministic case, the time dependence of the 
second-order response is determined by the converging 
series which contains the matrix elements A n (t; s). The 
expansion should be formally performed after the series 
summation. The long-time expansion of (^i , ^2) is 
an asymptotic, rather than a converging expansion in 
e~ 2tl and e -2 ' 2 . We have developed a method, equivalent 
to regrouping, which allowed to approximate the infinite 
sum by a sum of a finite number of terms. Due to the 
alternating character of the series, any smooth cut-off 
effective for larger term numbers did not influence the 
result. 

The noise k > actually introduces a smooth cutoff in 
the sum over n in Eq. (|71[) for a given pair of resonances 
A, fj,. The suppressive exponential factor is present for 
arbitrarily small positive k. The series over angular har- 
monics is almost alternating and converging. Therefore, 
the sum of the series does not depend on the value of k 
as long as it is small. It is intuitively clear that another 
form of regularizing noise would lead to the same results. 

For the second-order response function we have demon- 
strated the equivalence of the spectral decomposition in 
the limit of vanishing noise k — > with the asymptotic 
expansion in the case of deterministic dynamics. In short, 
this follows from the equality lim K _>o A n (t; s) — A n (t; s) 
for the evolution operator matrix elements and the prop- 
erty that the asymptotic expansion of S^^ii,^) origi- 
nates from the expansion of A n (t] s) and the subsequent 
proper summation of the apparently diverging series. 



VIII. DISCUSSION 

In the present manuscript we have studied linear and 
second-order nonlinear response of a stochastic system 
obtained by adding Langevin noise to a deterministic sys- 
tem whose dynamics is strongly chaotic. The determin- 
istic system is represented by a free particle moving on 
a 2D compact Riemann surface M 2 with constant neg- 
ative curvature (geodesic flow). We chose the random 
Langevin force to be orthogonal to the particle momen- 
tum, so that the noise does not change the energy, and 
random walk occurs on the energy shell represented by 
the 3D reduced phase space M 3 . The stochastic dynam- 
ics has been analyzed in terms of the Fokker-Planck oper- 
ator C(k) = —(k/2)& 2 + L, where the second term stands 
for the deterministic component (advection), whereas the 
first term describes the Langevin noise in the form of dif- 
fusion in the momentum space with k being the diffusion 
coefficient. 

The k = case that corresponds to purely deter- 
ministic dynamics has been studied in our earlier work 
0] , where we have employed strong dynamical symmetry 
(DS) in the system to find the analytical solution to the 
problem. The Langevin noise added to the system has 
been chosen in a way that it does not break down the 
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DS. Similar to the deterministic case, where the space of 
reduced phase space distributions has been decomposed 
into simpler components [irreducible representations of 
the DS group 50(2, 1)] invariant with respect to the Li- 
ouville operator L, the same components form invariant 
subspaces for our Fokker-Planck operator C{n). 

Here we have considered noise as a regularization 
to construct spectral decompositions of the Perron- 
Frobenius operator and focused on the noiseless K — ► 
limit of the stochastic dynamics. We have employed the 
DS to analyze the eigenvalue problem for the Fokker- 
Planck operator C whose eigenfunctions are nice smooth 
functions for finite values of k. In each irreducible com- 
ponent the distributions can be represented by functions 
on the circle, whereas the Fokker-Planck operator is rep- 
resented by a sum of the first-order Liouville operator 
and a ID diffusion operator — (n/2)d 2 . The determinis- 
tic dynamics has the stable u = 7r and unstable u = 
fixed points that represent the dynamical processes along 
the stable and unstable directions, respectively. Mapping 
the original problem onto much simpler II? stochastic 
dynamics on the circle allowed for an explicit analysis of 
the eigenvalue problem. Being interested in the noiseless 
limit we focused on small finite values n <C 1 and found 
the relevant eigenfunctions of C(n) analytically using the 
WKB method where the diffusion coefficient k plays the 
role of the square H 2 of the Planck constant. 

The fluctuation-dissipation theorem relates the linear 
response function to the two-point correlation function 
calculated earlier [ll[, the latter being interpreted as 
an expansion in RP resonances using the language of 
rigged Hilbert spaces. The long-time asymptotic expres- 
sion for the second-order response functions, has the form 
of a spectral decomposition over the same set of reso- 
nances [4]- Interpretation of the asymptotic expansion 
in the nonlinear case is more involved due to the follow- 
ing reasons. The eigenmodes that correspond to the RP 
resonances are represented by generalized, rather than 
smooth functions. In the linear case the initial smooth 
distribution should be decomposed in the RP modes. The 
signal is computed by convoluting the RP modes with the 
smooth dipole function. Both operations are well-defined 
for generalized functions. In the nonlinear case the sec- 
ond interaction with the driving field involves applying 
a differential operator to a generalized function followed 
by projecting it onto another generalized function. The 
legitimacy of the latter operation is less obvious, and has 
been related to the cancellation of apparently dangerous 
terms Q. 

Interpretation of nonlinear response in terms of RP res- 
onances has been addressed in this manuscript by consid- 
ering the noiseless limit of the Langevin dynamics. In the 
nonzero noise case k > 0, the spectral decomposition of 
the response functions of any order is conceptually more 
or less straightforward. We have demonstrated explic- 
itly for geodesic flows that in the n — > limit the rele- 
vant eigenvalues of C(k) converge to the RP resonances, 
whereas the coefficients in the spectral decompositions 



of the response functions converge to the coefficients of 
the asymptotic series in the purely deterministic k = 
expressions derived in Ref. 0. 

Summarizing, RP resonances can be interpreted as 
the noiseless limit for the eigenvalues of the Fokker- 
Planck operator C(k), and the spectral decompositions 
in this limit reproduce the long-time asymptotic series 
for the response functions. Note that the dynamical 
(■-function can be also reproduced as the limit C( z ) = 
lim re _>o Z^ 1 det(z — C(n)). In our model the RP decom- 
position for the linear response is represented by a con- 
verging series, whereas the nonlinear response is given 
by a non-converging asymptotic series. The converging 
character of the spectral decomposition is lost when the 
limit k — ► is applied. Computation of the expansion co- 
efficients in the nonlinear case requires a delicate summa- 
tion of almost sign-alternating series whose convergence 
is ensured thanks to the noise. 

Applications of our results are not limited to interpre- 
tations of purely deterministic quantities. Irreversibil- 
ity that shows itself in the decaying correlations appears 
only when the deterministic chaotic dynamics is regu- 
larized by some kind of coarse graining. Full "physical" 
mixing always requires some diffusion mechanism. The 
difference between stable and chaotic deterministic dy- 
namics is that the diffusion-induced "physical" mixing is 
much more efficient in a chaotic system. We consider a 
small fraction of the phase space that represent the initial 
conditions. As the ball of initial conditions is stretched 
and folded back, the shape becomes elongated along un- 
stable directions and contracted along stable ones, while 
the phase space volume remains constant. The diffusion 
time scales as a square of the blurring size. Therefore, a 
purely diffusive relaxation in the stable system occurs on 
the time scale of r rog ~ I 2 / n. In a chaotic system the scale 
of the density inhomogeneity decreases with time expo- 
nentially. Considering for simplicity a 3D reduced phase 
space of a chaotic Hamiltonian system, a ball of size a 
of initial conditions becomes a fettuccine-like shape. The 
fettucine length grows as ae xt . The size of the system I 
induces the folding of the unstable manifold. We estimate 
the number of folds as ~ ae Jl. Then the characteris- 
tic "physical" mixing time r c ~ A" 1 In (/ 3 a _2 (A/«;) 1 / 2 ) 
in a chaotic system with weak diffusion is very short 
compared to the stable case. Note that the "physical" 
mixing time r c can be measured in photon echo experi- 
ments, where it represents the characteristic time scale of 
the photon echo decay as a function of the delay between 
the exciting and the dynamics-reversing pulses. 



APPENDIX A: EIGENMODES OF THE 
FOKKER-PLANCK OPERATOR 

In this section we calculate the relevant eigenmodes 
of the Fokker-Planck operator in a given representation 
of 5*0(2, 1) labeled by s. They are given by symmetric 
solutions of the Schrodinger equation on a circle with the 
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effective Hamiltonian (|28 



k d 2 , . 
2du~ 2 ^ u) + 



sin 2 u 
2k 



scosm £\(u) = \£\{uXAl) 



The Hamiitonian is self-adjoint with respect to a natural 
symmetric scalar product. 

For our purposes we focus on the weak limit c < 1. 
The low lying states with |A| <C k^ 1 are concentrated 
near the the potential minima, in the vicinity of u = 
and u = tt. The principal approximation for the solutions 
of Eq. (|A1|) in the regions \u\ <C 1 or \u — tt\ <C 1 can 
be found by retaining up to second order terms in the 
expansion of the potential. 

For |u| < 1 we arrive at the approximate equation: 



and A\ and B\ are two complex constants. 

We begin with the construction of a symmetric solu- 
tion near u = it. The requirement of the solution to be 
invariant with respect to {u — tt) <-» [tt — u) combined 
with the relations between Hermite functions leads to an 
identity 



Bj_ 



T(-n)cos^ ' 



(A7) 



and after some transformations of the Hermite functions 
we arrive at an explicit form of a symmetric solution for 
(EH): 



£(u)= A 



2-^+1) e ^_ W2 x 



(A8) 



2 du 2 



£A(u) + ;r &(«) = (* + *)&(«)■ ( A2 ) 
2k 



that reproduces the well-known Schrodinger equation for 
a linear harmonic oscillator. One can easily identify the 
characteristic scales \/k/k = 1 and (kk) 1 / 4 = ^/k of 
the energy and length, respectively. A general solution 
of Eq. (|A2[) can be represented in terms of the Hermite 
functions as 



(A3) 

with Ao and £?o being complex constants. The parameter 
v is related to the energy A by 



X-l + s. 



(A4) 



As opposed to the case of a harmonic oscillator, the 
general solution (|A3|) contains both decaying and grow- 
ing waves. Their relative amplitude can be determined 
only by solving the equation on the whole circle. The 
functions H v (z) and e z H- v -\(iz) that are linearly in- 
dependent for any parameter v can be simply related 
to confluent hyperbolic functions and parabolic cylinder 
functions, the latter being the solutions of the original 
Schrodinger equation with the harmonic potential. We 
prefer to deal with the Hermite functions rather than 
with other special functions since the former reproduce 
Hermite polynomials in the case of integer order. 

The general solution of Eq. (|A2|) near u = tt has the 
form similar to Eq. (|A3[) : 



««) = A ie -(— r/^ Hui <^_^ 



Bie («- 7 r) 2 /2« iJ _ ;yi _ i 



(A5) 



where 



v\ = v — 2s = A — — — s , 



(A6) 



I U — 7T \ / U — TT 

H- Ul -i [ i — — ) + H- Vl -i I i 



Therefore, in the region \u — tt\ <C 1 the symmet- 
ric solution is determined by two unknown parameters 
v = v\ + 2s and A±; both can assume complex values. 

The harmonic approximation for the potential can be 
used if | it — tt\ <C 1, hence it is the region where the 
solution given by Eqs. (|A51 IA7|) is valid. 

The general WKB solution of Eq. (jAip under the bar- 
rier and around u = tt /2 is represented by a superposition 
of two waves: 



A 2 



exp(— S(u)) 



Bo 



exp(5(u)), (A9) 



where p(u) and S(u) are defined by 



j(u) = V sin 2 u — 2ks cos u — 2k\ , (A10) 



(All) 



S(u) — — dwp(w) 
« J 

it/2 



The quasiclassical expansion over k for the wavefunction 
phase can be employed when the wavelength k/p(u) does 
not change too fast: 



d 



du \p(u) 



<C 1 



(A12) 



This holds when we are not too close to the classical turn- 
ing points where p(u) turns to zero. Small value of k en- 
sures that the classical turning points lie close to u = or 
u = tt. The WKB solutions under the barriers are valid 
at least for \u\ 3> \/2k(\ + s) and \u — tt\ 3> ^/2k(A — s). 
The right-hand sides of these inequalities contain expres- 
sions for the classical turning points. If |A|, \s\ ~ 1, the 
validity of the quasiclassical approximation (|A9[) condi- 
tions is limited to the regions 



\u — tt\ > \[k . 



(A13) 



We can invoke a further approximation in the expres- 
sions for the solutions under the barrier. For k <C 1 and 
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| sinu| 3> \J~k~ (i.e. in the regions determined by the same 
inequalities (|A13| ). we arrive at a simplified expression 
for S(u), where only those terms that are small compared 
to one have been neglected: 



S(u) =- 



dw I sin w — ks 



cosu 



k; A 



7T/2 

COS It 



(A14) 
(A15) 



— s In sin u — A In tan — . 

2 



In approximating the pre-exponential factor of Eq. 
it is sufficient to set p(u) = sin it. 

We can see that in all intervals the wave function con- 
sists of two components that exponentially grow in the 
opposite directions under the barrier. A crude estimate 
for the magnitudes of the two components is presented 
in Fig. [T] Only the dominant part sin 2 uj (2k) of the po- 
tential contributes to the exponential dependence of the 
barrier transparency on k. 




FIG. 1: Schematic picture of the amplitudes that result from 
the tunneling through the potential barrier. 

We further notice that for it close to it (actually in the 
interval ft 1 / 2 <C (n — it) <C k 1 / 4 -C 1) we can use the 
expansion 



S(u) 



(U - 7r) 2 

2k 



+ (A-s)ln(7r-u)-Aln2. 



where the neglected terms are small compared to one. 

This is the interval where the WKB solution under 
the barrier and the solution in the harmonic well both 
represent a good approximation and can be matched. 
We use the asymptotic expressions for the Hermite func- 
tions valid for \u — ir\ 3> \[k. In fact, taking into ac- 
count Eqs. (|A4[ IA6[) . we can see that the term with 
A 2 in Eq. (|A9p and the solution given by Eq. (j3U|) con- 
tain the same quadratically increasing argument in the 
exponential, and the same power-like pre-exponential de- 
pendence: 



A 2 (tt -u)- x+s -h x exp 



1 

K 



(U - 7r) S 

2k 



(A16) 



Si 



2(tt - it) 



exp 



iit(vi + 1) (it — 7r) z 



2k 



In Eq. (I A16[) we have also used the symmetry £(it) 
£(27r — it), which implies that for positive (u — n) 3> 



the asymptotic form of £(it) 
that is proportional to Bi . 

Therefore, the matching near u 



36|) includes the only term 
7r yields 



B\ = A 



2<° 



(A17) 



The value of A\ is then found from Eq. (|A7|1 . The com- 
ponent of £(it) in Eq. (|A5|) with B 2 in front of it decreases 
rapidly cx exp(— cosit/t«) under the barrier as u deviates 
from 7T. Its magnitude can be estimated as 



B 2 cx A\e 



(A18) 



The procedure of solving the Schrodinger equation 
near it = is completely similar: In the region ^/H <S 
|u| <C 1 both approximate solutions (|A3[) and (|A9[) are 
adequate and can be matched. To that end we expand 
S'(ii) in the interval |u| <§; k 1 / 4 , neglecting the terms that 
are small compared to one: 

1 u 2 

S2(u) = h (A + s) In it + Aln2 , 

k 2k 

and also employ the asymptotic forms of the Hermite 
functions in the overlap region. With u approaching 
the second term in Eq. (|A9|) decreases and hence matches 
the second term in Eq. (|A3() . 



Bq cx B 2 e » cx A\e 



(A19) 



The component with A 2 in front of it matches under the 
barrier with the component that has Aq in front of it and 
decreases as u deviates from 0: 



-4, 



A e~ 



(A20) 



First we assume that the ratio B\j A\ given by Eq. (|A7[) 
does not contain exponentially small term ~ e -1 / K . Then 
Eqs. (|A17P - (|A20p imply that the second component in 
Eq. (|A3|) scales B cx e~ 4 / K A near it = 0. Its presence 
is related to the nonresonant tunneling along the circle 
from the potential minimum into itself. The argument in 
the suppressive exponential factor is evaluated as 



2tt 

— / dupiu) — 

K J K 




du I sinitl = 



(A21) 



The tunneling is nonresonant, since the two wells near 
u = and u = ir are offset by s which is not an integer 
number. 

Performing the same procedure in the region with 
sin it < as it has been done for sin it > 0, we find the 
symmetric solution 



Z(u)=A e-" 2 ^H„^ , 



(A22) 



valid for |ii| <^ 1. For arbitrary v the function may have 
a cusp at u = 0. The solution is smooth if v is a non- 
negative even integer. Then the Hermite functions H v 
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reduce to Hermite polynomials of even order. The de- 
pendence X u = v+ (1 — s)/2 on v with v = 0,2, ... of the 
eigenstate energy for the states concentrated near u = 
is given by Eq. (|A4[) . Exponentially small corrections 
oc e~ 4 / K , induced by the tunneling can be neglected. On 
the other hand, the power-like in n corrections to the 
energies, play a crucial role for v > 0. The latter correc- 
tions appear because of the deviation of trigonometric 
functions from their harmonic approximations. However 
the method developed in Section fVTl allows to avoid these 
calculations completely. 

Complex conjugation of the Hamiltonian that cor- 
responds to the change of sign of s, or to the shift 
u — > u + 7r implies the quantization condition v\ = 2k 
with k = 0, 1, 2, . . . for the other set of eigenstates. They 
have energies A = 2fc+(l+s)/2 and are concentrated near 
u = tt. The eigenfunctions satisfy the relations similar to 
the specified above. 

The eigenfunctions of the Hamiltonian ([2"8|) found 
above significantly differ from those of the linear har- 
monic oscillator only far from the interval |u| < k 1 / 4 . 
Outside the interval the exponentially decaying functions 
take negligibly small values so that the normalization fac- 
tor has the form typical to the harmonic oscillator: 

A 0}U = (^ K )- 1 / 4 (2^!)- 1 / 2 . (A23) 



The phase of the normalized eigenfunctions of Hermitian 
operators can be chosen arbitrarily. For our Hamiltonian 
and scalar product (|29|) this reduces to the freedom in 
the sign choice. 

Finally, for illustrative purposes, we calculate the cor- 
rection to the energy = v — s + 1/2 of the eigenstate 
that originates from the vicinity of u = ( note that the 
result is not used anywhere): 



— TT 

-(2s(l + 2i>)- (l + 2i/ + 2// 2 )). (A25) 



One of the ways derive Eq. (IA24I) is to consider the next 
term in the WKB expansion under the barrier. The 
WKB wavefunction contains the state energy. One can 
simply require that the wavefunction does not contain 
logarithmic terms due to their absence near u = 0. 
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